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We present a dynamical spectral model for Large Eddy Simulation of the incompressible magneto- 
hydrodynamic (MHD) equations based on the Eddy Damped Quasi Normal Markovian approxima- 
tion. This model extends classical spectral Large Eddy Simulations for the Navier-Stokes equations 
to incorporate general (non Kolmogorovian) spectra as well as eddy noise. We derive the model 
for MHD and show that introducing a new eddy- damping time for the dynamics of spectral tensors 
in the absence of equipartition between the velocity and magnetic fields leads to better agreement 
with direct numerical simulations, an important point for dynamo computations. 
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I. INTRODUCTION 

Magnetic fields permeate the universe. If kinetic effects 
such as Hall current, ambipolar drift of anisotropic pres- 
sure tensor, may be prevalent at small scales, the large- 
scales can be described in the magnetohydrodynamic 
(MHD) approximation. For example, electric fields and 
ionospheric currents play a dynamic role in the evolution 
of the atmosphere above 100 km, and the input of energy 
from the magnetosphere during magnetic storms can af- 
fect the thermosphere and ionosphere on global scales. 

MHD has many similarities with Navier-Stokes (NS) 
turbulence: recall the Batchelor analogy between vortic- 
ity and induction, both undergoing stretching through 
velocity gradients (see equation (j6j)). On that basis, one 
can conjecture that the energy spectrum will be of the 
Kolmogorov type, as in fact observed in numerical simu- 
lations of a decaying flow ([1] and references therein), as 
well as in the Solar Wind [2[. However, Iroshinkov and 
Kraichnan (IK hereafter) hypothesized that the slowing- 
down of nonlinear transfer by Alfven waves would alter 
the energy spectrum [3], predicting a spectrum ^ & -3 / 2 , 
as recently observed in numerical simulations 0, 0, 0] 
and in Solar Wind observations [7]. The Lagrangian 
renormalized approximation also gives spectra compat- 
ible with the IK model Q. These spectra are isotropic, 
but the presence of a strong quasi-uniform magnetic field 
B at large scale renders the dynamic anisotropic. One 
can compute exactly the reduced dynamics in that case 
0, using weak turbulence theory; the emerging energy 
spectrum ~ |k^| -2 , where refers to wave vectors per- 
pendicular to B (note that the isotropization of such 
a spectrum is compatible with the IK spectrum). Note 
also that a weak turbulence spectrum was observed in 
the magnetosphere of Jupiter [l0[ , the evidence stemming 



from an analysis of Galileo spacecraft data. A weak tur- 
bulence spectrum also obtained as well in a large numer- 
ical simulation of MHD turbulence in three dimensions 
|6[ at a magnetic Taylor Reynolds number of ~ 1700 at 
scales smaller than where the isotropic IK spectrum is 
observed. 

Both the terrestrial and Jovian magnetospheric plas- 
mas as well as the solar wind, the solar atmosphere, and 
the interstellar medium are highly turbulent conduct- 
ing compressible flows sustaining magneto-acoustic wave 
propagation and a better understanding of their dynam- 
ics, leading for example to star formation, requires ad- 
equate tools for modeling them. Furthermore, there is 
currently a surge of interest for achieving an experimen- 
tal dynamo (the growth of a seed magnetic field through 
fluid motions, see [ll|). In the case of liquid metals, or 
the fluid core of the Earth [l2|, EE 03 or the solar con- 
vection zone, the magnetic Prandtl number is very small 
(10 -5 or less); hence, the dynamo instability occurs in a 
turbulent flow and modeling the turbulence in order to 
study this phenomenon is in order [l5|, [l6[ • 

There are few models for MHD (see e.g. the recent 
review in [l7[), comparatively to the fluid case where 
the engineering community has been driving a vigorous 
research agenda. In that light, this paper aims at devel- 
oping such a model, in the context of a spectral approach 
following the work of Chollet and Lesieur [l8[ for the fluid 
case, using two-point closure of turbulent flows. We give 
the basic equations in the next section and then move 
on in Section IIIII to recall the EDQNM closure formu- 
lation for MHD. New triad relaxation times are intro- 
duced in Section [IV] and first tested in Section [VI The 
case of a random flow for two values of the magnetic 
Prandtl number Pm are treated respectively in Section 
ED^m = 1) and Section [VH] (P M = 0.1), and the deter- 
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ministic Orszag-Tang flow in three dimension is analyzed III. SPECTRAL MODELING 

in Section [Villi Finally, Section HXl is the conclusion and 

two technical appendices ([Aland lB]) are given at the end. j± The original EDQNM closure 



II. MAGNETOHYDRODYNAMIC EQUATIONS 

Let us consider the Fourier transform of the velocity 
v(x, t) and the magnetic B(x, t) at wavevector k: 



/oo 
v(x,t)e- ik - x dx 
-OO 

/OO 
B(x, t)e" ik - x dx. 
-OO 



(1) 
(2) 



The MHD equations describe the time evolution of a con- 
ducting fluid velocity field coupled to a magnetic field. 
They derive from Maxwell's equations with the assump- 
tion that velocities are sub-relativistic, hence the dis- 
placement current can be neglected [H, [2(3] • In terms of 
the Fourier coefficients of the velocity and the magnetic 
components, the MHD equations with constant unit den- 
sity read: 



d_ 

dt 

d_ 

dt 



uk 2 v(k,t) = t y (k,t) 



rjk 2 )b(k,r) = t M (k,i) , 



(3) 



(4) 



with k • v = in the incompressible case and k • b = 
indicating the lack of magnetic monopoles in the classical 
approximation. Here, b = B/y7^oPo is the Alfven veloc- 
ity, with fiQ the permeability, and po the uniform density 
(taken equal to unity); r] is the magnetic diffusivity, v the 
kinematic viscosity, and t v (k, t) and t M (k, i) are bilinear 
operators for energy transfer written as: 

t a v (k,t) = -zP a/3 (k)/c 7 ^2 v (3(Pi t ) v -/( c Li t ) 

p+q=k 

+iP a(3 (k)k 7 Yl hMb^t) (5) 

p+q=k 

t a M (k,t) = -iSapk^ ^2 bpiPit^fat) 

p+q=k 

M^Wp**) > ( 6 ) 

p+q=k 

with P a p(k) = S a /3 — k a kp/k 2 a projector that allows to 
take the pressure term of the velocity equation into ac- 
count via the Poisson formulation. The magnetic Prandtl 
number is defined as Pm = v/v- Finally note that, in 
the absence of dissipation [y — = 77), the total energy 
E T = 0.5 < v 2 + b 2 >, the correlation between the ve- 
locity and the magnetic field He =< v • b > and the 
magnetic helicity < A • b > (with b = V x A) are invari- 
ants of the ideal MHD equations in three dimensions. 



The Large Eddy Simulation model (LES) derived in 
(2lj (Paper I hereafter) is now extended to the MHD 
equations in its non- helical version (LES-P). As a first 
step, a spectral filtering of the equations is realized; 
this operation consists in the truncation of all veloctity 
and magnetic components at a wave-vector k such that 
|k| > k c where k c is a cut-off wave number. In an in- 
termediate zone lying between k c and 3k c both kinetic 
and magnetic energy spectra are assumed to behave as 
power-laws followed by an exponential decrease: 



E v (k,t) 



kp ^ k <C 3/C/- 



(7) 



E M (k,t) = E™k- a *e- d z\ k c <k<3k c , (8) 

where E^f , and a|f, #|f , E^f , are evaluated at 

each time step of the numerical simulations, through a 
mean square fit of the resolved kinetic and magnetic en- 
ergy spectra. In a second step one can write the modeled 
MHD equations as: 

[d t + (v(k\k c ,t) + vk 2 )]v a (k,t) = £<(k,t) (9) 
[d t +(r,(k\k c ,t) + V k 2 )]b a (k,t) = t%<(k,t) (10) 

where the < symbol indicates that the nonlinear transfers 
are integrated on the truncated domain such as p + q = k 
with |p|, |q| < k c . The quantities v (k\k Cl t) and r] (k\k Cl t) 
which are respectively called eddy viscosity and magnetic 
eddy diffusivity are expressed as: 



v{k\k c ,t)-- 



(SY(k,p,q,t) + SY(k,p,q,t)) A 
u ~t n t-~i t /" / 7 — \ -apaq 



r](k\k 



2k 2 E v (k,t) 



(S^(k,p,q,t) + Sf(k,p,q,t)) 
2k 2 E M (k,t) 



(11) 
dpdq 
(12) 



(see Paper I for more details). Here the S i ' (k,p,q,t) 
terms (see Appendix A), correspond to the absorption 
terms of the Eddy Damped Quasi Normal Markovian 
(EDQNM) nonlinear transfer, leading in particular to 
turbulent eddy diffusivities (see e.g. [22] for the MHD 
case). A > is the integration domain on /c, p, triangles 
such as p and or q are larger than k c and both p and q 
are smaller than 3k c . 

Finally, to take into account the effect of the emis- 
sion (eddy-noise) terms of the EDQNM nonlinear trans- 
fer (i.e. SY(k,p,q,t), SY(k,p,q,t), S^(k,p,q,t), and 
S3 4 (k,P,q,t)), we use a reconstruction field procedure 
which also enables to partialy rebuild the phase relation- 
ships between the three spectral components of each ve- 
locity and magnetic fields, as explained in Paper I [2lj . 
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B. First numerical tests 

We first implemented our LES model with the 
EDQNM equations of Pouquet et al [23] (see also Q). In 
this formulation the triad-relaxation time Skpq (see Ap- 
pendix A) takes three characteristic times into account: 
- a (combined) dissipation time td defined as: 



TB 1 (k) = (v + l 1 )k 2 

a nonlinear time rs expressed as: 



Ts\k) = A 







q 2 [E V (q)+E M ( q )]dq 



(13) 



(14) 



and an Alfven time ta which reads: 



(15) 



This constitutes a straightforward generalization of the 
EDQNM closure to the case of MHD flows (see for exam- 
ple |24j ) in which two new times, specifically the Alfven 
time and the diffusion time built on magnetic resistiv- 
ity, are incorporated in a phenomenological manner. A 
comparison of a simulation using this LES model (LES 
MHD I, or run II in Table D}, with a DNS simulation 
(run I in Table [Ij) is shown in Fig. [1] (see next section for 
more information on the numerical procedure). One sees 
that both the kinetic and magnetic energy spectra are 
overestimated by the model at scales close to the cut-off, 
indicative of an inadequate energy transfer in the model 
at these scales. When evaluating numerically the differ- 
ent eddy dampings using eqns. (|13H15p , one observes that 
the Alfven time is almost one order of magnitude shorter 
than all other times, including the diffusion time at the 
smallest resolved scales (see j25j); this leads to an insuf- 
ficient damping at the scales close to the cut-off. This is 
in part due to the dominance of the magnetic energy at 
large-scale; in that sense, it could be linked to the partic- 
ular flow under study and parametric analyses of several 
flows will have to be performed in the future in order 
to fine-tune this MHD model. However, the discrepancy 
displayed in Fig. [1] could also be linked to the particu- 
lar expression of eddy damping chosen in [23]. We are 
thus led to examine more closely the dynamics of energy 
transfer within the EDQNM framework. 



IV. NEW RELAXATION TIMES FOR EDQNM 

We now analyze the precise structure of the equations 
leading to the EDQNM closure; this is done in Appendix 
iBl Note that a similar but more complex and system- 
atic approach can be found in [8] in the context of the 
Lagrangian renormalization approximation. Our analy- 
sis results in the expression of three new eddy diffusiv- 




FIG. 1: Kinetic (top) and magnetic (bottom) energy spectra 
at time t = 1, £ = 3, £ = 5, and £ — 10 from upper to lower 
curves for data I (256 3 DNS, solid line), and data II (64 3 LES 
MHD I, plus symbols). 



ity times, with which one can build four different eddy- 
damping rates jikpq-, namely: 



,, vv - 

r'kpq 

VM 
l^kpq 



-T VV 
- T D 



\k,p,q) +r N l(k,p,q) , 



= Td M \k : p,q)+T N 1 L (k : p,q)+r A 1 (k : p,q) 



Z>^D V ~\^p, q )^r- 1 L (k,p, q )^rX 1 (Kp, q ) ,(16) 
\k,p,q) +T~l(k : p,q) + T^ 1 (fc,p,g) (17) 



..MM _ 

H'kpq ' 



MM 



where: 



T N l(k,p,q) = r N l(k)+T N l(p)+T N l(q) , (18) 



~}{k,p, q) 



- v \k)+T2\p)+T^\q) . (19) 



with t a (k^p^q) based on a new Alfven-like time, namely 



\j k EV(q)dq) 



[ k q 2 E M (q)dq 
Jo 



(20) 

Finally, from these eddy-damping rates we derive the 
corresponding triad-interaction times (as defined in Ap- 
pendix A). For the EDQNM kinetic energy equation, we 
obtain two different O's, ®Xpq a PP^ e( ^ t° SY \k,p,q,t) 
and SY(k,p,q,i), and Qf p f applied to S]f(k,p,q,i) and 
SY(k,p,q,t). 



For the EDQNM magnetic energy equation we also 
obtain two triad- interaction time, namely O^^f applied 
to Sf(k,p,q,t) and S^(k,p,q,t) and 6^ applied to 
S3 1 (k,p,q,t) and S^f (k : p,q,t). Note that this formula- 
tion leads to distinguish between the Joule and the vis- 
cous dissipation more systematically than in [23], and 
therefore to possibly better simulate flows at magnetic 
Prandtl numbers different from unity. Another differ- 
ence is the special treatment of the Lorentz force in the 
velocity equation. We call LES MHD II the model that 
takes these new triad-relaxation times into account. 



TABLE I: Parameters of the simulations. Initial conditions 
I.C., grid resolution A", kinematic viscosity z/, and magnetic 
Prandtl number Pm = v/r], with 77 the magnetic diffusivity. 
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DNS 


Random 


256 


2.e- 3 


1 


II 


LES I 


Random 


64 


2.e- 3 


1 


III 


LES II 


Random 


64 


2.e- 3 


1 


IV 


DNS 


Random 


256 


8.e- 4 


0.1 


V 


LES II 


Random 


64 


8.e- 4 


0.1 


VI 


LES I 


Random 


64 


8.e- 4 


0.1 


VII 


DNS 


OT 


256 


2.e- 3 


1 


VIII 


LES II 


OT 


64 


2.e- 3 


1 



V. NUMERICAL SET-UP 



In order to assess the ability of the model to reproduce 
the physics involved in MHD flows, we performed Direct 
Numerical Simulations (DNS) of the three-dimensional 
MHD equations, at a resolution of 256 3 grid points, to- 
gether with computations based on our LES MHD for- 
mulation, but now using 64 3 grid points. We performed 
this comparative study from three different simulations of 
freely decaying MHD flows. To test our model in a simple 
configuration, we first run a simulation at Pm = 1 with 
random initial conditions and no correlation between the 
velocity and the magnetic field (run I for the DNS, run 
II for the LES MHD I, and run III for the LES MHD 
II, in Table [J). Since the new eddy-damping times we 
derived allow for a clear distinction between the kine- 
matic viscosity and the magnetic dissipation, we then 
simulated a flow with identical initial conditions but with 
now Pm = 0.1 (run IV for the DNS, and run V for the 
LES MHD II in Table QJ. We recall that, in the work of 
[23^ the EDQNM closure has been derived for the case 
when the cross-correlation (or cross-helicity) spectrum 
H c (k) =< v(k) • b(k)(/c) > is assumed to be identi- 
cally zero in the presence of helicity (see p], [26| for the 
non helical case in the presence of velocity-magnetic field 
correlation). However, for many flows, this quantity is 
non negligible; furthermore, it can be strong locally (in 
particular in the vicinity of vorticity and current sheets) 
even when the global correlation is close to zero [2?], UK . 
We thus performed as well a simulation at Pm = 1 for 
which the velocity and the magnetic field are signifi- 
cantly correlated, in order to see how our model may 
adapt to such a situation. We chose the so-called three- 
dimensional Orszag-Tang flow (run VI for the DNS, and 
run VII for the LES MHD II in Table [TJ) for which ini- 
tially 2H C /E T = 0.5. 

From all these simulations, we studied global flow quan- 
tities such as the total, kinetic and magnetic energies, as 
well as helicities, and the cross-correlation energy. We 
also analyzed the spectral behaviors of these quantities. 



VI. RANDOM FLOW AT P M = 1 

We first investigate the model behavior for a flow 
with random initial conditions, presenting no cross- 
correlation, and at magnetic Prandtl number of unity. 



A. Inter-comparison of models 

In this section, we compare the efficiency between the 
model that involves the eddy damping times stemming 
from [23] (LES MHD I), and the model where the new 
eddy-damping times derived in Appendix B are included 
(LES MHD II). In Figure we plot the relative differ- 
ence of the kinetic and magnetic energy spectra com- 
puted from both LES models with the ones computed 
from the DNS. The spectra are chosen at time t = 1, 
close to the time of maximum dissipation. 

At large scales (bewteen k = and k ~ 15) the LES 
MHD II model globally gives a better approximation of 
the kinetic and magnetic energy spectra. Bewteen k = 14 
and k = 20 for the kinetic energy spectra, and between 
k = 15 and k = 25 for the magnetic energy spectra, the 
LES MHD I seems to give better results. This is due 
to the fact that LES MHD I and DNS spectra cross at 
a wavenumber located inside these ranges. Finally, at 
small scales, the LES MHD II data lead to a much better 
approximation than the LES MHD I data. At different 
times, the comparison between LES MHD I and LES 
MHD II results leads to similar results (not shown). We 
therefore focus our study on the LES MHD II model for 
the remainder of the paper. 



B. Global quantities 

We here study the time evolution of the global kinetic, 
E v (t), and magnetic, E M (t), energies for runs I and III, 
as shown in Figure [3j). 

One can observe that the modeled kinetic and mag- 
netic energies both closely follow the DNS evolutions, 
although at short times (between t = 1 and t = 5 for 
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FIG. 2: Lin- log plots of the relative difference with DNS en- 
ergy spectra for the velocity (top) and the magnetic field (bot- 
tom) at time £ = 1, for runs II and III, compared to the DNS, 
run I. Note the large error in LES I at large k. 
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FIG. 3: Temporal evolution of the kinetic and magnetic en- 
ergy for runs I (256 3 DNS), and III (64 3 LES MHD II). 



E M (t) and between £ = 1 and £ = 3 for E v (t)), the 
model slightly under-estimates them. 
Since our field-reconstruction procedure uses the flow (ki- 
netic and magnetic) helicities (even though the model 
itself does not take into account at this stage the heli- 
cal contributions to evaluate the transport coefficient), 
we plot in Fig. [4] the time evolution of both kinetic and 
magnetic helicities (respectively H v (t) and H M (£)). 

One can notice that, even though both modeled kinetic 
and magnetic helicities do not exactly match the DNS 
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FIG. 4: Time evolution of the kinetic and magnetic helicities 
for runs I (256 3 DNS), and III (64 3 LES MHD II). 



results at each time, they remain close and reproduce the 
main DNS time fluctuations. Note that the LES MHD I 
model provides similar results. 

We do not present here the temporal evolution of the 
cross-helicity H c '(£), since it is negligible when compared 
to the total magnetic and kinetic energy. Indeed, this 
correlation, initially equal to zero, reaches a maximum 
value of 0.081 for the DNS run, and of 0.069 for the LES 
MHD II run, to respectvely finish at a value of 0.051 and 
0.056. 

We now investigate the spectral behavior of our LES 
model by comparing the DNS and LES MHD II kinetic 
and magnetic energy spectra at various dynamical times. 




FIG. 5: Total energy spectra E T (k) = E M \k) + E v {k), at 
time £ = 1, £ = 3, £ = 5, and £ = 10 from top to bottom, for 
runs I (256 3 DNS solid line), and III (64 3 LES MHD II +). 

Figure [5] shows the total (kinetic plus magnetic) energy 
spectra E T (k) = E M (k) + E v (k) at times £ = 1, £ = 3, 
£ = 5, and £ = 10 obtained from DNS and LES MHD 
II computations. At any wavenumber and at any time, 
our LES MHD II model reproduces more correctly the 
DNS spectra than the LES MHD I does (see Fig. [TJ. It 
is clear that the spectral over-estimations at small scales 
obtained with this latter model is cured by the new for- 
mulation of the eddy-damping rates. 
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VII. RANDOM FLOW AT P M = 0.1 

io ( 



Since the new eddy- damping times involved in our LES 
MHD II model allow for a more refined differenciation 
between the magnetic diffusivity and the kinematic vis- 
cosity, we simulated a flow at a magnetic Prandtl number 
less than unity, namely Pm = 0.1. In order to highlight 
the efficiency of the new damping times to reproduce the 
flow dynamics, we compared both the LES MHD I and 
II data against the DNS results. For these simulations 
we kept identical flow initial conditions as in the previous 
section. A first comparison between the time evolution 
of the kinetic and magnetic energies computed from a 
DNS, and a simulation using the LES MHD II model, is 
plotted in Figure 




FIG. 7: Total energy spectra, at times £ = 1, £ = 3, £ = 5, 
and £ = 10 from top to bottom, for runs IV (256 3 DNS, solid 
line), V (64 3 LES MHD II, plusses), and VI (64 3 LES MHD 
I, triangles) at Pm =0.1. 
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FIG. 6: Total kinetic and magnetic energy temporal evolu- 
tion, for runs IV (256 3 DNS), and V (64 3 LES MHD II) at a 
magnetic Prandtl number of 0.1 



VIII. DETERMINISTIC ORSZAG-TANG FLOW 
AT P M = 1 

For a majority of flows, the correlation between the 
velocity and the magnetic field (or cross-helicity) is non 
negligible, leading to a slowed-down dynamics and energy 
spectra depending on the amount of correlation in the 
flow [29( . It has also been observed that local patches of 
either aligned or anti-aligned velocity-magnetic field con- 
figurations can be found both in the solar wind and in 
numerical simulations [2?], HH • We therefore decided to 
evaluate the ability of our model to treat a flow with 
strong cross-correlation by examining the evolution of 
the so-called three-dimensional Orszag-Tang flow with 
an initial global correlation H c (t = 0) = 1.63 (to be 
compared with the total kinetic and magnetic energy 
E v (t = 0) =E M (t = 0) = 2). 



One can observe that the model almost reproduces the 
exact temporal evolution of both kinetic and magnetic 
energy. The evolution of the kinetic and magnetic helic- 
ities (not shown) is also well-reproduced by the model. 
Once again, the cross-correlation remains weak all along 
the simulations; initially equal to zero, it reaches a maxi- 
mum value of 0.056 for the DNS, 0.057 for the LES MHD 
I, and 0.057 for the LES MHD II runs, before respective 
final values of 0.044 (DNS), 0.046 (LES MHD I), and 
0.045 (LES MHD II). 

We now present in Figure. [71 the total (kinetic plus 
magnetic) energy spectra evolution at times t = 1, t = 3, 
t = 5, and t = 10, obtained from DNS, LES MHD I, 
and LES MHD II data. Although at small wavenumbers 
both LES models correctly reproduce the DNS spectra, 
at large wavenumbers, strong differences appear among 
these various spectra. Indeed, the LES MHD II results 
slightly underestimate this range of the DNS spectra, 
whereas the LES MHD I highly overestimates it. 



A. Global quantities 



~ 1.5 





DNS:E v (t) 
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DNS: E M (t) 

LES MHD II: E M (t) _ 
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\ NX. 
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time 



FIG. 8: Kinetic and magnetic energy evolution, for runs VII 
(256 3 DNS solid line), and VIII (64 3 LES MHD II dashed 
line) with non-zero velocity-magnetic field correlation. 
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The kinetic energy evaluated with the LES MHD II 
fits with great accuracy to the kinetic energy obtained 
with the DNS (see Fig. [8]); however, the magnetic energy 
which is well-reproduced until £ = 2 departs measurably 
from the DNS data after this time. 




2 4 6 8 10 12 

time 



FIG. 9: Global velocity- magnetic field cross correlation, for 
runs VII (256 3 DNS straight line), and VIII (64 3 LES MHD 
II dashed line). 

The global cross-correlation, computed from either DNS 
or LES MHD II data, are quite close (see Fig. [9]), demon- 
strating that although the model does not explicitly 
take this quantity into account, it still maintains a re- 
liable evolution for it. However, the well-known tempo- 
ral growth of the normalized cross-correlation coefficient 
pit) = H c (t)/(E v (t) + E M (t)) shown in Fig. [TO] is not 
represented as accurately as either Et or H c . This could 
be tentatively attributed to the fact that turbulent trans- 
port coefficients based on the velocity-magnetic field cor- 
relation itself would emerge from a complete model (as 
derived in see also [26]) the effect of which might be to 
dampen the correlation growth over time. Note that this 
discrepancy likely emerges from the less accurate repre- 
sentation of the magnetic energy itself, as displayed in 

Fig. m 




time 



FIG. 10: Correlation coefficient p(£), for runs VII (256 3 DNS 
solid line), and VIII (64 3 LES MHD II dashed line). 



B. Spectral features 

We finally investigate the spectral behavior of our 
model on this particular Orszag-Tang flow. We respec- 
tively plot in Fig. [TP and Fig. [12] the kinetic and magnetic 
spectra of both DNS and LES MHD II, at times t = 1, 
£ = 3, £ = 5, and £ = 10. 




k 



FIG. 11: Kinetic energy spectra, at times £ = 1, £ = 3, £ = 5, 
and £ = 10 from up to down, for data VII (256 3 DNS straight 
line), and VIII (64 3 LES MHD II plusses). 
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FIG. 12: Same as Fig. [TT]for magnetic energy spectra. 

One can observe strong similarities of the modeled 
spectra with the directly simulated ones, although small 
differences appear at large scales. 

In order to evaluate the effect of the model on the cross- 
correlation, scale by scale, we represented in Fig. [13] its 
spectra at times £ = 3 and £ = 10 (we show only two 
times for readability purpose). 

We can observe that at large scales which are the most 
energetic, the model reproduces correctly the spectra ob- 
tained with the DNS at both times. But, close to the 
cut-off, the model strongly under-estimates the cross- 
correlation. This phenomenon, as stated before, is linked 
to the eddy viscosity and eddy diffusivity which dissipate 
the kinetic and magnetic resolved scales, as well as the 
cross-correlation at these scales. The reconstruction pro- 
cedure allows to reinject energy and helicity (when taken 
into account) at these scales, but not the correlation. 



8 



10" 



~ 10 

2C 



10" 





DNS 

+ LES MHD II 






+\ 


+ \ 



10 u 



10' 



FIG. 13: Correlation spectra at t = 3 (top) and t = 10 (bot- 
tom), for data VII (256 3 DNS straight line), and VIII (64 3 
LES MHD II +). 



IX. CONCLUSION 

In this paper we accomplish two complementary tasks. 
We first develop a LES for MHD using the EDQNM equa- 
tions and transport coefficients derived in [23[ but in the 
non- helical case. We then show that not all relevant time 
scales appearing in the cumulant expansions of the prim- 
itive MHD equations are taken into account in the phe- 
nomenological formulation of (23[. Indeed, one can derive 
several new eddy-damping times for the EDQNM equa- 
tions, and document how, by using them, one can consid- 
erably improve the treatment of the magnetic and kinetic 
energy transfers in the Large Eddy Simulation approach 
taken in this paper. 

A possible extension of this work is to be able to in- 
corporate the effect of either cross-helicity [26] and of ki- 
netic and magnetic helicity [23[ in the evaluation of eddy 
viscosities and eddy noise. The fact that the modeling 
algorithm does not depend on a specified inertial index 
may also be of some help in the case of a high velocity- 
magnetic field correlation when different spectra emerge 
at high values of the (normalized) H c cross-helicity [29j . 

Furthermore, with such a model many astrophysical 
and geophysical flows can be studied and perhaps more 
importantly a vast range of parameters, in particular the 
magnetic Prandtl number, can be examined. Among 
such problems, the generation of magnetic fields at ei- 
ther low or high magnetic Prandtl number is of prime 
importance, in particular in the former case in view of a 
set of laboratory experiments studying this effect [30[ • 



APPENDIX A: EDQNM CLOSURE 

For completeness, we recall here the expression of the 
EDQNM closure equations for the magnetic and kinetic 
energy without helicity. The first non-helical EDQNM 
equations were first derived in [22[ but we follow here 
the notation of [23[ which gives the free helical closure: 



(d t - 



2uk 2 )E v (k,t) = T v {k,t) 



2 V k 2 )E M (k,t) = T M (k,t) 



(Al) 
(A2) 



where the nonlinear transfer terms for the Jrinetic and 
magnetic energy, respectively T y (/c,t) and T M (k,t) are 
expressed as: 



'<*•*> - SL 



f M (k,t) 



kpil (t)S v (k,p,q,t)dpdq, (A3) 
kpq (t)S M (k,p,q,t)dpdq . (A4) 



Here is the integration domain with p and q such that 
(k,p,q) form a triangle, and 9 k (t) namely the triad- 
relaxation time is expressed as: 

Ok P q(t) = , (A5) 

f^kpq 

with Hkpq = fJ>k + + where the /i&'s are called eddy- 
damping rates and read: 

+ y|fc^V^ +(v + v)k 2 . (A6) 

The constant A can be expressed as a function of the 
Kolmogorov constant Ck appearing in front of the kinetic 
energy spectrum such that: 

A = 0.218Cf , (A7) 

following [l8[. 

The expressions of S v (k,p,q,t) and S M (k,p,q,t) can 
be further explicited (with the time dependency of mag- 
netic and kinetic energy spectra omitted here) as : 

S V (k,p,q,t) = ^b kpq [k 2 E V {q)E V {p)-p 2 E V (q)E V {k)} 

+ j q c kpq [k 2 E M (q)E M (p)-p 2 E M (q)E V (k)] 

= SX{k,p,q,t)+S%(k,p,q,t) 

+ SX(k,p,q,t)+SX(k,p,q,t) . (A8) 
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S M {k^q,t) = ^h kpq [k 2 E M (p)E v (q)-p 2 E v (q)E M (k)} 

-c kpq \%E v {p)E M (q)-E M {q)E M (k) 
pq lp z 

SF(k,p,q,t)+S?{k,p,q,t) 



fc 3 



S^(k,p,q,t)+S^(k,p,q,t) . 



(A9) 
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In Eqs. (|A8|) and (|A9ft the geometric coefficients bk pq , 
Ckpq, and hkpq are defined as: 

h Pq = pk~ x (xy + z 3 ), c kpq = pk' 1 z(l - y 2 ) , 
h kpq = z(l - y 2 ) , 

where x,y, z are the cosine of the interior angles opposite 
to k, p, q. This completes the description of the EDQNM 
closure for MHD as developed in [22|, [23[. The helical 
case, dealt with in [2l[ for a pure fluid and in [23[ from 
the EDQNM standpoint, will be studied in a forthcoming 
paper when coupling to a magnetic field is involved. 



APPENDIX B: 



A MORE GENERAL EDDY 
DAMPING 



In [23[, the eddy damping term is built on a phe- 
nomenological ground; namely, one argued about the ne- 
cessity of introducing the Alfven time scale in the damp- 
ing coefficient, without actually referring to the set of 
cumulant expansion equations. This change alone, from 
a traditional hydrodynamic EDQNM closure, led to en- 
ergy spectra that differ from the Kolmogorov case, with 
a k~ z l 2 law in the uncorrelated case, and with a E ± (k) ~ 
k~ m in the correlated case [1], with m + +m~ = 3; here, 
E ± (k) are the energy spectra of the Elsasser variables 
z ± = v±b. 

However, when examining the succession of equations 
for the higher-order moments, and keeping the total cor- 
relation between the velocity and the magnetic field equal 
to zero to simplify the algebra, a more complex structure 
emerges, which may help the modeling of the MHD dy- 
namics to be closer to the DNS than the results shown in 
FigHJ There are in fact four groups of terms in the clo- 
sure equations. The first group corresponds to the pure 
fluid case and can be written symbolically as: 



+ v(k 2 + p 2 + q 2 ) I < uuu >~ k < uuuu > . 



d_ 
dt 

(Bl) 

This leads, as usual, to two characteristic times written 
here as: r^ v = (v(k 2 + p 2 + q 2 )) 1 and tnl = (ku) -1 . 
The second group writes symbolically again as: 



dt 



+ vk + r](p +<?)]< uuu > 



(B2) 



k < uuuu > +fc < bbuu > . 



Here, two new times can be extracted, namely a dissi- 
pative time r^ M = {vk 2 + rj(p 2 + q 2 )) , and ta = 
u(kbb) _1 ^ a modified Alfven time. 

The third group of closure terms is of the following 
type: 

d \ 
— + (rjk 2 + rjp 2 + vq 2 )\ < bbu > (B3) 



~k< bbbb > +k < bbuu > . 
and finally the fourth group: 
d 



dt 



vp 



r\q l ) <bbu> 



(B4) 



k < bbbb > +k < bbuu > . 



Again, the following new characteristic dissipative 
times can be a priori deduced from these two groups 
(with similar nomenclatures as before): t^ m = 



( V k 2 



r\p* + vq" 



and 



r MV 



( V k 2 



vp" 



In conclusion, a careful examination of the cumu- 
lant equations has led to the adjunction of several new 
times, distinguishing between magnetic and kinetic en- 
ergy transfer as well as the different quantities entering 
the transfer terms. 

Note that the new modified Alfven time is finally ex- 
pressed as: 



[\k)=C A 



J k E^(q)dq y\r 
J k EV(q)dq) [Jo 



q 2 E M {q)dq 



(B5) 

It incorporates the lack of equipartition between the ki- 
netic and magnetic energy that is often observed, and 
this for example should also alter the dynamics, in the 
early (kinematic) phase of the dynamo problem. The 
non-linear time has the classical expression built only on 
the velocity field: 



NL 



q 2 E v 



(q)dq 



(B6) 



Finally, we numerically estimated the value of the Alfven 
time constant Ca = 0.8. This point will need further 
study as we extend the number of flows that are tested 
with this LES. The constant A is determined through 
the relation |A7j in Appendix A. The model has thus 
two open parameters that can be evaluated once the con- 
stants appearing in front of the energy spectra in MHD 
are determined. 

Also note that the way the dissipation coefficients 
are taken into account may well affect the results when 
the magnetic Prandtl number differs substantially from 
unity unless possibly when both the kinetic and magnetic 
Reynolds numbers are very large because of the effect of 
renormalisation of transport coefficients [3l| . 

Similarly to equation (|A5|) for the eddy damping rate 
in [23[, we define generalized rates as: 



1 



-u XY t 



u XY 



(B7) 



fi^ Y and with XY standing 



with fijg 

for either VV, VM, MV or MM and with: 

= (T% v (k)y 1 + (T NL (k)y 



VV 



(B8) 
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d M = {r v D M (k))- 1 + (r NL (k))- 1 + (r A (k))- 1 , (B9) 



tf v = {rg v {k)) 1 + (r NL (k))- 1 + (7 A (k))- 1 , (BIO) 
and 

r = (t% M (k))- 1 + (r NL (k))- 1 + (r A (k))- 1 . 

(Bll) 
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